El análisis predictivo emplea datos históricos para predecir eventos futuros. Normalmente, los datos históricos se utilizan para crear un modelo que capture diversos patrones temporales. Este modelo predictivo se usa entonces con los datos actuales para predecir lo que acontecera. Existen diversas metodologias para llevar acabo predicciones, en este taller optamos por el enfoque bayesiano para llevar acabo la prediccion.
El enfoque bayesiano, en el contexto del analisis predictivo, enfoca su analisis en el calculo de la distribucion predictiva aposteriori. Esta distribucion emplea la informacion de los datos asi como la distribucion aposteriori de los parametros de interes . Empleando estos insumos la distribucion predictiva podemos generar datos futuros; es decir, datos faltantes en un determinando horizonte.
\[\begin{align*} f(y_{pred}|y_{obs}) = \int_{\theta} f(y_{pred},\theta|y_{obs})d\theta =\int_{\theta} f(y_{pred}| \theta y_{obs})f(\theta|y_{obs})d\theta \end{align*}\]
Actualmente existen diversas metodologias para llevar acabo el calculo numerico de esta distribucion las cuales en R corresponden al empleo de un determinado paquete en. Algunos paquetes ampliamente conocidos son los siguientes :
Relativamente hace no mucho tiempo la metolodogia de la Aproximación anidada integrada de Laplace o INLA(por sus siglas en ingles) ha cobrado importancia por su relativa eficiencia para estimar modelos complejos.
El analisis de prediccion como en la practica puede ser entendida como predecir datos desconocidos en un determinado horizonte. En ese sentido estos datos a futuro pueden ser interpretados como datos faltantes o missings values, los cuales son el objeto de la predicción.
| Año | Semana | Numero de muertes |
|---|---|---|
| 2019 | 1 | 0 |
| 2019 | 2 | 1 |
| 2019 | 3 | 4 |
| 2020 | 1 | 0 |
| 2020 | 2 | 1 |
| 2020 | 3 | 4 |
| 2021 | 1 | NA |
| 2021 | 2 | NA |
| 2021 | 3 | NA |
Para llevar acabo predicciones empleando INLA, una aproximacion es justamente la anterior. Una vez dispuestos los datos como en la tabla anterior, basta con ajustar un modelo tipico en INLA, ya que INLA internamente calculara la distribucion predictiva aposteriori .
library(tidyverse)
library(INLA)
library(yardstick)
library(gt)
library(innovar)
library(spdep)
library(reshape2)
db <- readRDS("db_excess_proc.rds")
db.lima <- db %>%
filter(prov=="LIMA")
db.lima %>%
ggplot()+
geom_line(aes(x=date,y=n),color="#011f4b")+
geom_point(aes(x=date,y=n),color="black")+
facet_wrap(vars(sexo))+
ylab("Numero de muertes")+
xlab("")+
ggtitle("Evolucion del Numero de muertes por género")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5,size = "22"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) Para llevar acabo el ajuste de modelos empleamos la funcion inla(), la cual tiene diversos parametros, algunos de los mas importantes son :
data : Un objeto tipicamente de la clase data.frame el cual representa los datos para ajustar el modeloformula : Un objeto de la clase formula que especifica la ecuacion que pretendemos ajustar como por ejemplo inla(target ~ 1 + x1). En la formula podemos especificar efectos lineales (introduciendo el nombre de la variable) o no lineales (empleando f()).verbose : Una variable del tipo boolean, la cual señala si se desea mostrar en consola el proceso de convergenciadb.femenino <- db %>%
filter(sexo=="FEMENINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))
db.masculino <- db %>%
filter(sexo=="MASCULINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))
db.n.femenino <- db.lima %>%
filter(sexo=="FEMENINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))
db.n.masculino <- db.lima %>%
filter(sexo=="MASCULINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))lima.m.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.masculino
)
lima.f.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.femenino
)Los parámetros antes detallados son los esenciales para ejecutar el ajuste de un modelo empleando INLA.Sin embargo, algunos parametros extra a tener en consideracion son los siguientes :
family:Objeto de clase character.Este parametro es crucial, pues determina la distribucion de la variable objetivo, por defecto se encuentra en family=Gaussian.control.compute:Objeto de la clase list permite especificar el calculo de criterios de informacion tales como aic,dic,waic.lima.m.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.masculino,
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
)
lima.f.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.femenino,
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
)En el analisis de prediccion (forecast) tipicamente se emplea datos en el formato de series de tiempo para llevarla acabo. Para ello comunmente se asume que las series de tiempo pueden ser descompuestas del siguiente modo :
\[\begin{align*} y_{t} = S_{t}+T_{t}+e_{t} \end{align*}\]
donde :
\(S_{t}\) : Estacionalidad
\(T_{t}\) : Tendencia
\(e_{t}\) : Ruido blanco
Los componentes antes detallados pueden ser modelados en el contexto de INLA empleando distribuciones apriori para cada componente tales como los siguientes:
variablef(variable, model = "ar1")f(variable, model = "rw1")f(variable, model = "rw2")Dado el comportamiento que exibe la series que prentedemos modelar, para los componentes :
#inla.list.models("likelihood")
lima.m.m1 <- inla(n ~ 1 + f(annus,model = "ar1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
verbose = F,
family = "nbinomial",
data = db.n.masculino
)
lima.m.m2 <- inla(n ~ 1 + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.masculino
)
lima.m.m3 <- inla(n ~ 1 + f(annus,model = "ar1") + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.masculino
)fit.m.m1 <- lima.m.m1$summary.fitted.values$mean
fit.m.m2 <- lima.m.m2$summary.fitted.values$mean
fit.m.m3 <- lima.m.m3$summary.fitted.values$mean
n.m <- db.n.masculino$n
datos <- list("AR(1)"=fit.m.m1,"RW(1)"=fit.m.m2,"AR(1)RW(1)"=fit.m.m3) %>%
as.data.frame() %>%
melt( variable.name = "modelo",
value.name = "fit") %>%
group_by(modelo) %>%
mutate(actual = n.m,
date = db.n.masculino$date) %>%
ungroup() %>%
mutate(min_inter = c(lima.m.m1$summary.fitted.values$`0.025quant`,
lima.m.m2$summary.fitted.values$`0.025quant`,
lima.m.m3$summary.fitted.values$`0.025quant`),
max_inter = c(lima.m.m1$summary.fitted.values$`0.975quant`,
lima.m.m2$summary.fitted.values$`0.975quant`,
lima.m.m3$summary.fitted.values$`0.975quant`)
)
datos %>%
ggplot() +
geom_line(aes(x=date,y=actual)) +
geom_line(aes(x=date,y=fit),color="#011f4b") +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter),
alpha=0.4, #transparency
fill="#011f4b") +
facet_wrap(vars(modelo)) +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted",
color = "black", size=0.7)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) #inla.list.models("likelihood")
lima.f.m1 <- inla(n ~ 1 + f(annus,model = "ar1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.femenino
)
lima.f.m2 <- inla(n ~ 1 + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
verbose = F,
family = "nbinomial",
data = db.n.femenino
)
lima.f.m3 <- inla(n ~ 1 + f(annus,model = "ar1") + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.femenino
)fit.f.m1 <- lima.f.m1$summary.fitted.values$mean
fit.f.m2 <- lima.f.m2$summary.fitted.values$mean
fit.f.m3 <- lima.f.m3$summary.fitted.values$mean
n.f <- db.n.femenino$n
datos <- list("AR(1)"=fit.f.m1,"RW(1)"=fit.f.m2,"AR(1)RW(1)"=fit.f.m3) %>%
as.data.frame() %>%
melt( variable.name = "modelo",
value.name = "fit") %>%
group_by(modelo) %>%
mutate(actual = n.m,
date = db.n.femenino$date) %>%
ungroup() %>%
mutate(min_inter = c(lima.f.m1$summary.fitted.values$`0.025quant`,
lima.f.m2$summary.fitted.values$`0.025quant`,
lima.f.m3$summary.fitted.values$`0.025quant`),
max_inter = c(lima.f.m1$summary.fitted.values$`0.975quant`,
lima.f.m2$summary.fitted.values$`0.975quant`,
lima.f.m3$summary.fitted.values$`0.975quant`)
)
datos %>%
ggplot() +
geom_line(aes(x=date,y=actual)) +
geom_line(aes(x=date,y=fit),color="#011f4b") +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter),
alpha=0.4, #transparency
fill="#011f4b") +
facet_wrap(vars(modelo)) +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted",
color = "black", size=0.7)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.m <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.m,
estimate = fit)
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Hombres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Hombres | |||
| modelo | mae | rmse | msd |
| AR.1. | 40.66904 | 50.09989 | 27.70859 |
| RW.1. | 34.64344 | 42.68089 | 27.66907 |
| AR.1.RW.1. | 30.38785 | 37.24076 | 27.67888 |
perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.f <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.f,
estimate = fit)
tbl.yrd.f %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Mujeres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Mujeres | |||
| modelo | mae | rmse | msd |
| AR.1. | 30.25786 | 35.90613 | 0.10282570 |
| RW.1. | 24.32470 | 30.10790 | 0.06330288 |
| AR.1.RW.1. | 18.04854 | 22.13154 | 0.07311218 |
El enfoque espacio-temporal, implica complementar el enfoque anterior complejizando el modelo tomando en cuenta la dimesion espacial. Para ello requerimos incluir efectos aleatorios espaciales para controlar por tal dimension de analisis. Cabe mencionar que modelos espaciales areales empleando INLA, requieren necesariamente la creacion de un grafo(graph) o una estructura de datos espaciales que captura la disposicion de las areas como una matriz de pesos espaciales.
\[\begin{align*} y_{t} = S_{t}+T_{t}+\nu_{t}+u_{t}+e_{t} \end{align*}\]
donde:
\(\nu_{t}\) : efectos estructurados
\(u_{t}\) : efectos no estructurados
Los 2 nuevos componentes \(\nu_{t},u_{t}\) pretenden tomar en consideracion en la prediccion la correlacion espacial presente en los datos . Al respecto :
Los efectos no estructurados corresponden a efectos aleatorios que pretenden controlar caracteristicas no observadas de cada area bajo estudio. Para modelarlas podemos emplear f(area, model = "iid")
Los efectos estructurados son efectos aleatorios que toman en cuenta explicitamente la estructura espacial
f(area, model = "besag")f(area, model = "besagproper")peru.m.m5 <- inla(n ~ 1 + f(annus,model = "linear") +
f(week,model = "rw1") +
f(id.sp, model = "besag", graph = W.peru),
data = db.m.sp,
control.compute = list(dic = TRUE,
waic = TRUE,
cpo = TRUE),
control.predictor = list( link = 1)
)peru.f.m5 <- inla(n ~ 1 + f(annus,model = "ar1") +
f(week,model = "rw1") +
f(id.sp, model = "besag", graph = W.peru),
data = db.f.sp,
control.compute = list(dic = TRUE,
waic = TRUE,
cpo = TRUE),
control.predictor = list( link = 1)
)perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.m <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.m,
estimate = fit)
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Hombres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Hombres | |||
| modelo | mae | rmse | msd |
| AR.1. | 30.25786 | 35.90613 | 0.10282570 |
| RW.1. | 24.32470 | 30.10790 | 0.06330288 |
| AR.1.RW.1. | 18.04854 | 22.13154 | 0.07311218 |
perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.m <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.m,
estimate = fit)
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Mujeres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Mujeres | |||
| modelo | mae | rmse | msd |
| AR.1. | 30.25786 | 35.90613 | 0.10282570 |
| RW.1. | 24.32470 | 30.10790 | 0.06330288 |
| AR.1.RW.1. | 18.04854 | 22.13154 | 0.07311218 |